# coding: utf-8

import logging
from math import log, pi, sqrt

try:
    from scipy.integrate import quad
    from scipy.special import k0
    from scipy.stats import norm
    
    NORM_DIST = norm(0.0, 1.0)
except ImportError:
    pass

MAX_BESSEL_VALUE = 500.0

def AbsBessel2Norm(abs_x):
    if abs_x <= 4:
        integr_value, err = quad(k0, 0, abs_x)
        prob = 0.5 - integr_value / pi
    else:
        if abs_x > MAX_BESSEL_VALUE:
            logging.warn('Bessel2Norm: abs(x) = %s', abs_x)
            abs_x = MAX_BESSEL_VALUE
        integr_value, err = quad(k0, abs_x, MAX_BESSEL_VALUE*2)
        prob = 2 / pi * integr_value
    assert err < 0.001
    return NORM_DIST.isf(prob)

def SignedSafeCall(func, x):
    if x != 0:
        res = func(abs(x))
        if x < 0 and res > 0:
            return -res
        else:
            return res
    else:
        return 0.0

def generate_Bessel2Norm():
    d = {}
    for n in xrange(1, int(MAX_BESSEL_VALUE)):
        precision = 10**(2 - int(log(n)/log(10)))
        assert precision > 0
        for p in xrange(0, precision):
            x = n - p / float(precision)
            print x
            d[x] = AbsBessel2Norm(x)
    return d

Bessel2Norm = lambda x: SignedSafeCall(AbsBessel2Norm, x)

# значения Bessel2Norm от 0.1 до 100.0 c шагом 0.1
BesselNormValues = [0.000, 0.276, 0.451, 0.593, 0.715, 0.824, 0.923, 1.015, 1.100, 1.180, 1.256, 1.329, 1.397, 1.464, 1.527, 1.588, 1.648, 1.705, 1.761, 1.815, 1.868, 1.919, 1.969, 2.018, 2.066, 2.112, 2.158, 2.203, 2.247, 2.291, 2.333, 2.375, 2.416, 2.457, 2.496, 2.536, 2.574, 2.612, 2.650, 2.687, 2.723, 2.525, 2.563, 2.601, 2.638, 2.675, 2.711, 2.747, 2.783, 2.818, 2.852, 2.887, 2.920, 2.954, 2.987, 3.020, 3.052, 3.084, 3.116, 3.148, 3.179, 3.210, 3.240, 3.271, 3.301, 3.330, 3.360, 3.389, 3.418, 3.447, 3.476, 3.504, 3.532, 3.560, 3.588, 3.615, 3.643, 3.670, 3.697, 3.724, 3.750, 3.776, 3.803, 3.829, 3.854, 3.880, 3.906, 3.931, 3.956, 3.981, 4.006, 4.031, 4.055, 4.080, 4.104, 4.128, 4.152, 4.176, 4.200, 4.224, 4.247, 4.270, 4.294, 4.317, 4.340, 4.363, 4.385, 4.408, 4.430, 4.453, 4.475, 4.497, 4.519, 4.541, 4.563, 4.585, 4.607, 4.628, 4.650, 4.671, 4.692, 4.714, 4.735, 4.756, 4.777, 4.797, 4.818, 4.839, 4.859, 4.880, 4.900, 4.921, 4.941, 4.961, 4.981, 5.001, 5.021, 5.041, 5.060, 5.080, 5.100, 5.119, 5.139, 5.158, 5.177, 5.196, 5.216, 5.235, 5.254, 5.273, 5.292, 5.310, 5.329, 5.348, 5.366, 5.385, 5.403, 5.422, 5.440, 5.459, 5.477, 5.495, 5.513, 5.531, 5.549, 5.567, 5.585, 5.603, 5.621, 5.638, 5.656, 5.674, 5.691, 5.709, 5.726, 5.744, 5.761, 5.779, 5.796, 5.813, 5.830, 5.847, 5.864, 5.881, 5.898, 5.915, 5.932, 5.949, 5.966, 5.982, 5.999, 6.016, 6.032, 6.049, 6.065, 6.082, 6.098, 6.114, 6.131, 6.147, 6.163, 6.179, 6.195, 6.212, 6.228, 6.244, 6.260, 6.276, 6.291, 6.307, 6.323, 6.339, 6.355, 6.370, 6.386, 6.402, 6.417, 6.433, 6.448, 6.464, 6.479, 6.494, 6.510, 6.525, 6.540, 6.556, 6.571, 6.586, 6.601, 6.616, 6.631, 6.646, 6.661, 6.676, 6.691, 6.706, 6.721, 6.736, 6.751, 6.766, 6.780, 6.795, 6.810, 6.824, 6.839, 6.854, 6.868, 6.883, 6.897, 6.912, 6.926, 6.940, 6.955, 6.969, 6.983, 6.998, 7.012, 7.026, 7.040, 7.055, 7.069, 7.083, 7.097, 7.111, 7.125, 7.139, 7.153, 7.167, 7.181, 7.195, 7.209, 7.223, 7.236, 7.250, 7.264, 7.278, 7.291, 7.305, 7.319, 7.332, 7.346, 7.360, 7.373, 7.387, 7.400, 7.414, 7.427, 7.441, 7.454, 7.467, 7.481, 7.494, 7.507, 7.521, 7.534, 7.547, 7.560, 7.574, 7.587, 7.600, 7.613, 7.626, 7.639, 7.652, 7.665, 7.678, 7.691, 7.704, 7.717, 7.730, 7.743, 7.756, 7.769, 7.782, 7.795, 7.808, 7.820, 7.833, 7.846, 7.859, 7.871, 7.884, 7.897, 7.909, 7.922, 7.934, 7.947, 7.960, 7.972, 7.985, 7.997, 8.010, 8.022, 8.035, 8.047, 8.059, 8.072, 8.084, 8.097, 8.109, 8.121, 8.133, 8.146, 8.158, 8.170, 8.182, 8.195, 8.207, 8.219, 8.231, 8.243, 8.255, 8.268, 8.280, 8.292, 8.304, 8.316, 8.328, 8.340, 8.352, 8.364, 8.376, 8.388, 8.399, 8.411, 8.423, 8.435, 8.447, 8.459, 8.471, 8.482, 8.494, 8.506, 8.518, 8.529, 8.541, 8.553, 8.564, 8.576, 8.588, 8.599, 8.611, 8.623, 8.634, 8.646, 8.657, 8.669, 8.680, 8.692, 8.703, 8.715, 8.726, 8.738, 8.749, 8.761, 8.772, 8.783, 8.795, 8.806, 8.817, 8.829, 8.840, 8.851, 8.863, 8.874, 8.885, 8.896, 8.908, 8.919, 8.930, 8.941, 8.952, 8.964, 8.975, 8.986, 8.997, 9.008, 9.019, 9.030, 9.041, 9.052, 9.063, 9.074, 9.085, 9.096, 9.107, 9.118, 9.129, 9.140, 9.151, 9.162, 9.173, 9.184, 9.195, 9.206, 9.216, 9.227, 9.238, 9.249, 9.260, 9.270, 9.281, 9.292, 9.303, 9.314, 9.324, 9.335, 9.346, 9.356, 9.367, 9.378, 9.388, 9.399, 9.410, 9.420, 9.431, 9.441, 9.452, 9.463, 9.473, 9.484, 9.494, 9.505, 9.515, 9.526, 9.536, 9.547, 9.557, 9.568, 9.578, 9.588, 9.599, 9.609, 9.620, 9.630, 9.640, 9.651, 9.661, 9.672, 9.682, 9.692, 9.702, 9.713, 9.723, 9.733, 9.744, 9.754, 9.764, 9.774, 9.785, 9.795, 9.805, 9.815, 9.825, 9.835, 9.846, 9.856, 9.866, 9.876, 9.886, 9.896, 9.906, 9.916, 9.927, 9.937, 9.947, 9.957, 9.967, 9.977, 9.987, 9.997, 10.007, 10.017, 10.027, 10.037, 10.047, 10.057, 10.067, 10.076, 10.086, 10.096, 10.106, 10.116, 10.126, 10.136, 10.146, 10.156, 10.165, 10.175, 10.185, 10.195, 10.205, 10.214, 10.224, 10.234, 10.244, 10.253, 10.263, 10.273, 10.283, 10.292, 10.302, 10.312, 10.321, 10.331, 10.341, 10.351, 10.360, 10.370, 10.379, 10.389, 10.399, 10.408, 10.418, 10.427, 10.437, 10.447, 10.456, 10.466, 10.475, 10.485, 10.494, 10.504, 10.513, 10.523, 10.532, 10.542, 10.551, 10.561, 10.570, 10.580, 10.589, 10.599, 10.608, 10.617, 10.627, 10.636, 10.646, 10.655, 10.664, 10.674, 10.683, 10.693, 10.702, 10.711, 10.721, 10.730, 10.739, 10.748, 10.758, 10.767, 10.776, 10.786, 10.795, 10.804, 10.813, 10.823, 10.832, 10.841, 10.850, 10.860, 10.869, 10.878, 10.887, 10.896, 10.905, 10.915, 10.924, 10.933, 10.942, 10.951, 10.960, 10.969, 10.979, 10.988, 10.997, 11.006, 11.015, 11.024, 11.033, 11.042, 11.051, 11.060, 11.069, 11.078, 11.087, 11.096, 11.105, 11.114, 11.123, 11.132, 11.141, 11.150, 11.159, 11.168, 11.177, 11.186, 11.195, 11.204, 11.213, 11.222, 11.231, 11.240, 11.248, 11.257, 11.266, 11.275, 11.284, 11.293, 11.302, 11.310, 11.319, 11.328, 11.337, 11.346, 11.355, 11.363, 11.372, 11.381, 11.390, 11.399, 11.407, 11.416, 11.425, 11.434, 11.442, 11.451, 11.460, 11.468, 11.477, 11.486, 11.495, 11.503, 11.512, 11.521, 11.529, 11.538, 11.547, 11.555, 11.564, 11.573, 11.581, 11.590, 11.598, 11.607, 11.616, 11.624, 11.633, 11.642, 11.650, 11.659, 11.667, 11.676, 11.684, 11.693, 11.701, 11.710, 11.719, 11.727, 11.736, 11.744, 11.753, 11.761, 11.770, 11.778, 11.787, 11.795, 11.804, 11.812, 11.820, 11.829, 11.837, 11.846, 11.854, 11.863, 11.871, 11.880, 11.888, 11.896, 11.905, 11.913, 11.922, 11.930, 11.938, 11.947, 11.955, 11.963, 11.972, 11.980, 11.988, 11.997, 12.005, 12.013, 12.022, 12.030, 12.038, 12.047, 12.055, 12.063, 12.072, 12.080, 12.088, 12.096, 12.105, 12.113, 12.121, 12.129, 12.138, 12.146, 12.154, 12.162, 12.171, 12.179, 12.187, 12.195, 12.203, 12.212, 12.220, 12.228, 12.236, 12.244, 12.252, 12.261, 12.269, 12.277, 12.285, 12.293, 12.301, 12.309, 12.317, 12.326, 12.334, 12.342, 12.350, 12.358, 12.366, 12.374, 12.382, 12.390, 12.398, 12.406, 12.415, 12.423, 12.431, 12.439, 12.447, 12.455, 12.463, 12.471, 12.479, 12.487, 12.495, 12.503, 12.511, 12.519, 12.527, 12.535, 12.543, 12.551, 12.559, 12.567, 12.575, 12.583, 12.590, 12.598, 12.606, 12.614, 12.622, 12.630, 12.638, 12.646, 12.654, 12.662, 12.670, 12.677, 12.685, 12.693, 12.701, 12.709, 12.717, 12.725, 12.733, 12.740, 12.748, 12.756, 12.764, 12.772, 12.780, 12.787, 12.795, 12.803, 12.811, 12.819, 12.826, 12.834, 12.842, 12.850, 12.858, 12.865, 12.873, 12.881, 12.889, 12.896, 12.904, 12.912, 12.920, 12.927, 12.935, 12.943, 12.951, 12.958, 12.966, 12.974, 12.981, 12.989, 12.997, 13.005, 13.012, 13.020, 13.028, 13.035, 13.043, 13.051, 13.058, 13.066, 13.074, 13.081, 13.089, 13.096, 13.104, 13.112, 13.119, 13.127, 13.135, 13.142, 13.150, 13.157, 13.165, 13.173, 13.180, 13.188, 13.195, 13.203, 13.210, 13.218, 13.226, 13.233, 13.241, 13.248, 13.256, 13.263, 13.271, 13.278, 13.286, 13.293, 13.301, 13.308, 13.316, 13.324, 13.331, 13.339, 13.346, 13.353, 13.361, 13.368, 13.376, 13.383, 13.391, 13.398, 13.406, 13.413, 13.421, 13.428, 13.436, 13.443, 13.450, 13.458, 13.465, 13.473, 13.480, 13.488, 13.495, 13.502, 13.510, 13.517, 13.525, 13.532, 13.539, 13.547, 13.554, 13.562, 13.569, 13.576, 13.584, 13.591, 13.598, 13.606, 13.613, 13.620, 13.628, 13.635, 13.642, 13.650, 13.657, 13.664, 13.672, 13.679, 13.686, 13.694, 13.701, 13.708, 13.715, 13.723, 13.730, 13.737, 13.745, 13.752, 13.759, 13.766, 13.774, 13.781, 13.788, 13.795, 13.803, 13.810, 13.817, 13.824, 13.832, 13.839, 13.846, 13.853, 13.860, 13.868, 13.875, 13.882, 13.889, 13.897, 13.904, 13.911, 13.918, 13.925, 13.932, 13.940, 13.947, 13.954, 13.961, 13.968, 13.975, 13.983, 13.990, 13.997, 14.004, 14.011, 14.018, 14.025, 14.033, 14.040, 14.047, 14.054, 14.061, 14.068]

def Bessel2NormFix(x):
    ai = min(1000, abs(int(x*10)));
    if (x >= 0):
        return BesselNormValues[ai]
    else:
        return -BesselNormValues[ai]
